

NASA CONTRACTOR REPORT 166307 


NASA-CR-1 66307 

OOID020 

Conservative Supra-Characteristics Method 
For Splitting the Hyperbolic Systems o£ 

Gasdynamics for Real and Perfect Gases 


For Preference 




NOT TO BE TA«Ef< moW ;>!lS ROOM 


C.K. Lombard 


LIBRARY BBPY 

IJ'R 1 fi 1982 


LANGLEY RESEARCH CENTER 
LIBRARY, NASA 
HAMPTON, VIRGINIA 


CONTRACT NAS2-10744 
January 1982 


niASA 


NF02311 




NASA CONTRACTOR REPORT 166307 


Conservative Supra-Characteris tics Method 
For Splitting the Hyperbolic Systems of 
Gasdynamics for Real and Perfect Gases 


C.K. Lombard 
PEDA Corporation 
1150 Fife Street 
Palo Alto, CA 94301 


Prepared for 

Ames Research Center 

under Contract NAS2-10744 


rvi/vsA 

National Aeronautics and 
Space Administration 

Ames Research Center 

Moffotl Field. California 94035 








Abstract 


A new conservative flux difference splitting is presented for the 
hyperbolic systems of gasdynamics. The stable robust method is suitable 
for wide application in a variety of schemes, explicit or implicit, 
iterative or direct, for marching in either time or space. The splitting 
is modeled on the local quasi one dimensional characteristics system 
for multi-dimensional flow similar to Chakravarthy 's nonconservative 
split coefficient matrix method (SCM); but, as the result of maintaining 
global conservation, the method is able to capture sharp shocks correctly. 
The embedded characteristics formulation is cast in a primitive variable 
the volumetric internal energy (rather than the pressure) that is effec- 
tive for treating real as well as perfect gases. Finally the relationship 
of the splitting to characteristics boundary conditions is discussed and 
the associated conservative matrix formulation for a computed blown wall 
boundary condition is developed as an example. 

The theoretical development employs and extends the notion of Roe 
of constructing stable upwind difference formulae by "sending" split 
simple one sided flux difference pieces to appropriate mesh sites. The 
developments are also believed to have the potential for aiding in the 
analysis of both existing and new conservative difference schemes. 





Acknowledgements 


The author wishes to acknowledge many stimulating discussions of 
Roe's method and upwind difference schemes with Professor Joseph Oliger 
of Stanford University who served as a consultant. Equally helpful 
has been the enthusiastic collaboration of Jaw Yen Yang in checking 
derivations and in programming and testing the method in model prob- 
lems. During the course of this research, to which he is making addi- 
tional contributions for his dissertation, Yang has been a Stanford 
Ph.D. candidate under Professor Daniel Bershader. Many earlier dis- 
cussions of the effectiveness of characteristic boundary conditions 
and Lax-Wendroff type hybrid characteristic-conservative methods with 
Thomas Coakley of NASA-Ames Research Center were a strong motivating 
influence behind the work reported here. 





RESEARCH REPORT 


CONSERVATIVE SUPRA-CHARACTERISTICS 
METHOD FOR SPLITTING THE HYPERBOLIC 
SYSTEMS OF GASDYNAMICS 
FOR REAL AND PERFECT GASES 

by 

C.K. Lombard 


I. Introduction 

For many years, Morretti^’^’^ (see also Pandolfi and Zanetti^) has 

advocated and demonstrated the physical content and accuracy available 

in the nonconservative primitive equations of gasdynamics and particularly 

the characteristics formulation of those equations, especially at bound- 

5 6 

aries. Following closely upon the Morretti " X scheme" Chakravarthy ’ , 
^ came forth with a simplified multidimensional, quasi one dimensional 
method based on eigenvector splittings of the coefficient matrices of 
the nonconservative Euler equations. Chakravarthy emphasized the cor- 
respondence between stable upwind differencings of the split equation 
pieces and the transformed quasi one dimensional characteristic 
relations. Getting generally very attractive results including very sharp 
captured shock transitions, Chakravarthy encountered only one substantial 
problem, that of unreliable shock location. 

On the other hand, a theorem of Lax and Wendroff^ that a finite dif- 
ference scheme in conservation law form can realize the Rankine-Hugoniot 
jump conditions automatically has led to a long line of development fea- 
turing such popular schemes as those of MacCormack^’^, Beam and Warming^^, 
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and Briley and McDonald^^. However, such conservative shock capturing 
schemes have been found deficient in resolution, tending to smear the 
transitions over several mesh points, and are subject to stability 
plaguing overshoots and undershoots and nonphysical solutions all of 
which must be controlled by the addition of ad hoc artificial dissipa- 
tion terms. Finally the artificial dissipation coupled with indiscrimi- 
nant directional differencing has led to excessive smearing of contact 
discontinuities. 

Recent recognition of these drawbacks has led to the introduction 

of directional signal propagation information in conservative schemes. 

The principal lines of development appear to be flux splitting intro- 

12 13 

duced first be Steger and Warming ’ and more recently flux difference 
splitting, put forth by Enquist and Osher^^, Osher and Solomon^^, and 
by Roe^^. Most recently Reklis and Thomas^^ have presented an upwind 
control volume centered method that appears to be a hybrid of flux split- 
ting and flux difference splitting. With the exception of the latter, 

18 

these methods have been surveyed and appraised by van Leer and by Marten, 
Lax and van Leer . Flux splitting suffers from the flaw that it fails 
to provide consistent directional signal information where it is most 
needed, namely in the vicinity of a change in eigenvalue sign. Flux dif- 
ference splitting on the other hand admits a greater generality than flux 
splitting and can be taylored to a variety of signal propagating philos- 
ophies or requirements as partially expressed in the works to date and the 
present report. 

The approach taken in the present work is to unify in so far as pos- 
sible the virtues of the signal propagation of the characteristics method 
with simple conservative upwind finite difference schemes in a way that 


2 



permits both explicit and implicit numerical methods. Among motivations 
for choosing the characteristics based splitting presented here is the 
belief that it is the natural and most numerically robust linearization 
of the gasdynamic equations since it supports the physical directional 
propagation of sound waves. This belief also seems to be supported by 
various gedanken experiments performed on numerical flow simulation pathol- 
ogies conceived out of past experience with conservative methods. 

In the characteristics interpretation and approach, the present method 
is closest to Chakravarthy; and in the development of general conservative 
upwind schemes out of one sided differences, follows Roe. Explicitly, 
reliance is placed on Roe's property U which the new matrix shares with 
his, though they are very different both in detail and concept of construc- 
tion. The principal feature of property U for admitting Rankine-Hugoniot 
satisfying shock transitions is 


aF = A Aq (1) 

where aF expresses the flux difference on either side of the transi- 
tion and Aq is the corresponding jump in the conservative variable 
vector. We find it more than mildly interesting that conservative methods 
based on this matrix, namely, 

6q + A Aq = 0 (2) 

are in the same form for method of solution as in their linear stability 
analysis. 
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II. Conservative Supra-Characteri sties Method 

Our method for constructing A rests on the observation that dif- 
ferences of nonlinear quantities such as pu appearing in the conservative 
equations can be expressed exactly as 

A pu = 7 AU + IT Ap . (3) 

This fact was recognized by Roe but the tougher triad terms evidently 
caused him to take a different ad hoc approach which he noted was 
general izable to other hyperbolic systems. The key to sorting out how 
to construct a splittable A from simple end point averagings ( ) of 
quantities across an interval is to analyze in detail how the charac- 
teristic equations are formed of the nonconservative primitive equations, 
on the one hand, and how the conservative equations are formed out of 
the primitive ones, on the other. The results will be sketched below 
but the cumbersome details are left to the Appendix. Here, we only 

remark that A is definitely not A , which van Leer has tried and 

19 

regards as too dissipative . 

Conservation And The Elemental Interval Equation 
Before describing the present A further, we introduce some relevant 
background material. The first item is the general discrete quasi one 
dimensional simple interval difference equation as an element of a multi- 
dimensional globally conservative numerical procedure. This concept is 
regarded as central to Roe's approach but is not introduced by him in 
this way. In the text figure below we see a bounded one dimensional 
computational domain divided up by interior mesh points into N-1 simple 
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intervals. We place data on the computed boundaries which are treated 
naturally with characteristic boundary conditions by this method, as will 
be described in a later section. 

f (1) • (2) . , , , . (N-1) \ 

123 N-1 N 

The simple (j)th interval equation whose splittings form the building 
blocks of the nodal point difference equations is then 

(J 6^q)j + (A^F)j. = 0 . (4) 

Here J is the volume of a slab from a computational topological cube 
bounded by all the adjacent mesh points to one at the center and isolated 
by passing a ^ constant coordinate surface through the center point. 

In one dimension J is just ax. . In multidimensions the "cube" (area 

J 

in 2-D) is bisected into slabs for each of the general curvilinear coor- 
dinate directions. For each slab and coordinate direction there exist 

•k 

elemental equations of the form (4) In the equation the term 6 q. 

T J 

is a measure of the change in the conservative variable q vector over 
the interval a? in the slab as a result of the bounding inviscid 
flux difference a F . (Note the discrete computational flux F 
contains the associated slab face area . For a discussion of 2-D 

and 3-D discrete finite volume formulations of the conservative equations 
in generalized curvilinear coordinates on a finite difference mesh see 

po pi 

Lombard, Davy, and Green and Thomas and Lombard . ) When such equa- 
tions (4) are summed over all mesh intervals in all coordinate directions. 
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all the interior flux terms cancel identically , and there results the 
discrete approximation to the global conservation law for the computational 
domain 


Jy q dV 3 ^ J 6^5 = -I F^. n s F.n dS . (5) 

Thus, when difference equations are constructed among split pieces of the 
elemental difference equations (4), global conservation requires only 
that the weighting functions for the distribution of each piece among 
the nodal difference equations should sum to one for each elemental equa- 
tion. 

While Roe focused his attention exclusively on just the flux term 
of equation (4), it is true that all discrete quasi-one-dimensional con- 
servative numerical methods can be constructed from the elemental dif- 
ference equations (4). Within this context, the difference between methods 

resides not just in the treatment of the flux term v^F but also in 

★ 

the definition of the unsteady term 6^q . For example for all methods 

•k 

constructed of simple one sided differences, q can be represented as 

^^j ^ “*^j ■ “^^j+i • 

In MacCormack's method, for instance, a = i . We will see presently 
that the way methods treat characteristic information is, or ought to be, 
reflected in the weighting function a . 

Stable Scalar Convection Equations 
The second item is the scalar convection equation which is well 
known to admit stable explicit and implicit numerical schemes when 
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backward or forward differenced according as the sign of the eigenvalue 
is positive or negative as below 

6^ Uj + A 7^ Uj = 0 , X ^ 0 

(7) 

6^ Uj. + X Uj = 0 , X £ 0 

Supra-Characteristics Matrix Splitting 
To continue with the discussion of the matrix that linearizes the 
flux difference equation (1), we are able to express A as 

A = M T D T"^ M"^ (8) 

Here the matrices M , T and T“^ bear considerable formal corre- 

12 13 

spondence to the matrices of flux splitting * but M is very 
different from the inverse M"^ . Specifically, is the matrix 

operator that constructs the right hand side (spatial) difference vector 
of the nonconservative primitive equations from the associated con- 
servative variable vector differences. The matrix T"^ is the operator 
that transforms the primitive equations into the scalar characteristic 
equations. Indeed, with reference to equation (4), we can write the 
1-D elemental characteristics equations that underlie the present 
conservative method as the set 

J D T'^ M'^ 6 q + D T"^ A,q = 0 (9) 

T 5 

In these equations, the diagonal matrix D is not the eigenvalue 
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matrix but a truth function which for a complete set of eigenvalues is 
the identity matrix. This formulation permits treating all the differ- 
ence terms of the primitive equations exactly without source terms and 
has led us to conceive the nomer conservative "supra(above, akin to)- 
characteristics" method (CSCM). To effect stable upwind difference 
methods the elemental characteristics equations are divided up for the 
jth interval according to partitioning the truth function D for positive 
and negative eigenvalues into and D" with complementary ones and 
zeroes on the diagonals. In analogy with equations (7) the stable one 
sided first order splittings of equations (9) are 

(J T’" M-")j 6^q.^j + (O'" T'" M’" A^q)^ = 0 

(10) 

(J D" T"^ M'^). 6 q. + (O' T"^ A o) . = 0 

J T J 5 J 

Note for the elemental interval j the forward difference A q. is 

s o 

equal to the backward difference • 

The stable conservative elemental difference equation pieces for 
the interval j are derived from equations (10) by multiplying by the 
matrix product M T . These, respectively from the right, serve to 
transform characteristics equations to primitive equations and primitive 
equations to conservative equations. The resulting splitting of equation 
(4) can be written 


(0 A^)j (A" = 0 

(J A-)j S^q. + (A- A^q)j = 0 


( 11 ) 
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Comparing these with equations (4) makes it evident that for a character- 
istics based splitting the elemental interval measure of the conservative 
variable vector q. for the unsteady term is defined (see (6)) as 

J 

q^ = A qj + A q^.^j . (12) 

With equations (11) as the basic stable elements and in analogy with 
Roe's prescription for "sending" the flux difference pieces A'*’ Aq and 
A" Aq , a host of stable conservative upwind methods can be constructed. 

In particular, a stable first order difference equation based on 
the elemental interval splitting (11) is 

[(J + (J A’)j] + (^ A^q)j_i + (A’ A^q)^ = 0 (13) 

Equation (13) is amenable to both explicit numerical methods involving 
a block diagonal matrix inversion and fully implicit methods involving 
block tri diagonal matrix inversions. Similarly a stable second order 
difference equation is 

i[-(J A‘^)j_2 + 3(J A'^)j._j + 3(J I'). - (J ^)j+i3 6^qj. 

+ K-(A^ A^q)j_2 + 3(A^ ^C^^j+1 ~ ^ 

Equation (14) is also amenable to explicit numerical methods and to 
implicit methods involving either block pentadi agonal inversions or block 
tridiagonal inversions, if the implicit convective terms are treated 
first order as in equation (13). Appropriate explicit and implicit 
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numerical methods and their stability for eigenvector split schemes 
have been discussed in references 12 and 13. 

For the elemental difference intervals adjacent to boundaries the 
distribution of the splittings (11) of the elemental interval equations 
to form stable interior nodal point difference equations such as (13) 
leaves pieces 


(J A")^6^qj + (A‘ A^q)j = 0 

(J A )fj_j (A 0 


(15) 


unused. These elemental equation pieces which embody the outrunning 
characteristics to the left and right boundaries respectively can be 
used as natural elements to construct dissipative boundary point dif- 
ference schemes that satisfy, with our interior point schemes, suf- 
ficiency conditions for stability of the coupled discrete initial boundary 

22 

value problem, Oliger . It is physically satisfying that the inclusion 
of these residual elemental equation pieces (15) for boundary point 
difference approximations is also necessary to close the global con- 
servation property of the set of elemental difference equations to the 
computational boundaries. 

The use of various akin space-time extrapolations from the interior 
has been discussed and analyzed for stability with different interior 

oo 24 

point schemes by Oliger and Sundstrbm and Gustafson and Oliger . 

25 

Yee , in a report with a recent extensive survey, sketched an approach 
based on flux splitting that is similar to the method we now present. 
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III. Natural Characteristic Boundary Conditions 

The diagonal matrices for interior nodal point difference equations 

such as (13) and (14) have a complete set of independent eigenvectors 

corresponding to the complete set of stable characteristics equations 

potentially available from either the right or left of the solution 

point depending on eigenvalue sign. At a computational boundary nodal 

point, however, only the stable subset of characteristics differenced 

from the interior and with negative eigenvalues for a left boundary 

and positive eigenvalues for a right, are applicable. Then, if there 

are N dependent variables in the conservative variable vector and 

only M stable characteristic relations available from the interior, 

we must specify N-M auxiliary boundary condition relations among the 

dependent variables. Such numerical auxiliary boundary conditions reflect 

in number and kind the natural boundary conditions for the PDE. To extend 

26 

the approach of Kentzer to a conservative matrix method such as we have 
here, auxiliary boundary relations are time differenced and included in 
the NxN block matrix of the diagonal to replace the unavailable exte- 
rior characteristic relations for a left boundary and A" for a 
right boundary. Without the introduction of these auxiliary boundary 
condition relations the split matrices A~ are singular. 

The method is illustrated for an isothermal wall boundary with pre- 
scribed blowing. Experience has shown it is most natural to introduce 
the auxiliary boundary conditions at the level at which the characteristic 
equations are formed namely at the T~^ matrix level. We specialize the 
construction to 2-D with a lower wall boundary which we assume is an n 
constant curvilinear coordinate line. 
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From the Appendix, the matrix T"’’ for the full set of characteristic 
relations is 


1 

0 

0 

_ 1 _ 

F 



w 

0 

5 

s 

0 


pc 

PC 






0 



X 


^c. 

PC 

YP 

0 





pc 

pc 



In the formulation this matrix is to be multiplied by the primitive equa- 
tion time difference (column) vector ( 6p , p 6u , p 6v , 6P . 

(Here for generality and applicability to real gas flows P is not the 
pressure p but the related .) The resulting characteristic vari- 
able differences are by row: (1) the entropy, (2) the tangential velocity, 
and (3) and (4) are the P"*^ and P“ normal pressure velocity compati- 
bility relations respectively. The eigenvalues associated with these 
characteristics are W’^ , W’^ , + nc , - ric . For subsonic blow- 

ing (inflow) only the latter of these is negative and is associated with 
a stable difference relation from the interior. Thus, the first three 
characteristic difference relations associated with the first three rows 
of T'^ need to be replaced by auxiliary boundary conditions. Such 
boundary conditions do not involve differences to the interior but only 
relations which serve to constrain the mutual variations of the dependent 
variables at the boundary. 

For consistency and global conservation with the interior point 
procedure, the boundary point computational procedure having the properties 
just described can conveniently be represented as 
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J M T(D"T'^M"^ + D‘^T'M')jSq^ + (M T D“T"^M’^A^q)^ = 0 . (16) 

Here for simplicity, first order differencing of the P' characteristic 
has been assumed from the wall. The auxiliary boundary conditions are, 
quite evidently, embodied in the matrix where is the truth 

function complement to D" , t' is the matrix in which theboundary 
conditions are formulated, and permits the formulation to be flex- 
ibly constructed among either the primitive variable differences, the 
conservative variable differences, or more commonly a combination of the 
two. 

Not too surprisingly it has been found that appropriate (stable) 
auxiliary boundary conditions have T' matrix representations row by 
row quite similar to the matrix rows of T~^ they naturally replace. 

Indeed for good matrix conditioning in the coupled difference equation 
boundary condition procedure, it is highly desireable that the corre- 
spondence be as close as possible. Consequently, in the following we 
elucidate the properties of T' in terms of the rows of T"^ that 
are naturally replaced. 

The first row of T~^ (multiplied by the primitive variable vector 
differences) can be seen to represent a logrithmic differential relation 
of a thermodynamic function the entropy. Indeed, a viable auxiliary^ 
boundary condition that can be applied is constant entropy (6 ln(_^) = 0), 
i.e. no wall heating, in which case row one of T"^ could be used intact 
for J' . A more general thermodynamic relation at the wall is the poly- 
trope law P ~ p“ which can be logarithmically differenced with the matrix 
coefficients 


0 


0 
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_ 1 _ 

ctP • 



For a = Y this is the isentropic wall and for a = 1 this is the 
isothermal or constant temperature wall. 

The second row of T"^ , for the tangential velocity equation, 

involves the direction cosines for the local tangent to the n constant 

coordinate lines. Thus at the wall the condition that blown flux enter 

the flow normal to the wall can be expressed as the inner product (the 

tangential velocity component) n' u - n' v = 0 . When differentiated 

y ^ 

this relation yields the same matrix representation for J' as the 
second row of T“^ . The third row of T’^ , for the pressure 
normal velocity compatibility equation, contains in columns two and 
three the direction cosines n' and nj for the local normals to the 
n constant coordinate lines. Many ablators of current and potential 
use sublimate at a nearly constant temperature, and the mass flow rate 
is proportional to the heating rate with a coefficient weakly dependent 
on the pressure. For convenience the pressure dependence can be described 
locally logrithmioally. Then the auxiliary boundary condition for the 
mass flux can be written as a variation about a local reference heating 
state as 


5 - Tp 



(17) 


Here m = pu + n" pv = pw"^ and the left hand side takes the T 
X y 

matrix representation 



p^n gp 
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Row one of T' is written for the primitive variable differences 
6p and 6P , while rows two and three are for the conservative variable 
differences 6(pu) and 6(pv) . The fourth row of T' is immaterial 
since it is not selected by the truth function and in fact need not 
be programmed at all. The M" that constructs the above desired com- 
bination of primitive and conservative variable differences from the homo- 
geneous conservative variable difference vector can be chosen as 

1 0 0 0 

0 10 0 

0 1 0 

u -V 1 


0 

u^ + v^ 
2 


IV. Results 

An indication of the potential effectiveness of the CSCM method can 
be found in the results of two model problems. 

The first model problem is the scalar Burgers equation ramp wave 
with inflow boundary condition and computed outflow. Here the single 
eigenvalue corresponding to A of a system is simply U and the equa- 
tion solved is 


6u^. + u^._^ v^.u = 0 for u > 0 . (18) 

In Figure la. the computed solution (boxes) is shown at three times 
(0.67, 1.33, 2.0) on the interval .5 _< x < 1 for initial conditions 
u = X on the interval. The exact solution u = plotted in 
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triangles connected by solid line has been given by Yee, Beam, and 
27 

Warming with computations that can be compared Figure lb. The solu- 
tion obtained with the present method was run explicitly at a Courant 
number of 0.8 and can be seen to be effectively exact. In particular 
there is no evidence of disconnectedness on the expansion indicating 
the method possesses a desirable degree of natural numerical dissipation. 

The second model problem is Shubins hyperbolic tangent nozzle prob- 
lem with supersonic inflow and subsonic outflow. The author is indebted 
to Warming and Yee for assistance in this problem by providing Shubins 
data and a code in which the method for the quasi 1-D Euler equations 

could conveniently be incorporated. Solution (Figures 2b. and 2c.), for 

25 

this problem with different methods have been published by Yee and 

27 

by Yee, Beam, and Warming . In Figure 2a. a solution is shown at 1000 
steps from an explicit calculation made with the present method at a 
Courant number of 0.8. For practical purposes the solution was obtained 
in 400-500 steps. The results, shown in boxes, can be seen to correspond 
to Shubins effectively exact solution to within the resolution of his 
data, shown connected by solid line. Again no difficulty is found with 
the expansion but also no noticeable smoothing of the jump transition is 
apparent. 

A comment on accuracy, the graphs have been labeled "first order 
upwind" after the finite difference convention of simple one sided dif- 
ferences employed for these results. However, owing to the fact that 
the method's eigenvalues are effectively centeree (the averaging proce- 
dure), the CSCM method is second order from a method of characteristics 
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point of view, Rakich . In any case, within the remarkable framework 
of Roe, the present finite difference method is extendable to any order 
of finte difference accuracy based on Taylor series. 

V. Summary 

A new conservative quasi one dimensional flux difference splitting 
has been presented for the hyperbolic terms of multidimensional gas- 
dynamics. The splitting based on local linearizing combinations of non- 
conservative primitive difference equations for one curvilinear coordinate 
direction is akin to a 1-D unsteady method of characteristics. In this 
respect the method is philosophically in accord with the split coefficient 
matrix method of Chakravarthy. But, unlike Chakravarthy and in accord 
with Roe, emphasis is placed on developing a method that maintains the 
global conservation property available in methods based on simple one sided 
conservative flux differences. To this end and extending the theoretical 
foundations of Roe's method, an elemental simple interval difference equa- 
tion is defined out of which any globally conservative finite difference 
method can be constructed. Akin to how the conservative PDE's are con- 
structed from the nonconservative primitive equations, a matrix transform 
of a discrete set of primitive difference equations is constructed to 
realize identically the conservative elemental interval equations. For 
application to real gas as well as perfect gas flows, the volumetric 
internal energy is chosen as the basis of the primitive equations repre- 
sentation rather than the pressure. 


★ 

private conversation 
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Neither the primitive difference equations nor their transform, the 
characteristic equations, are solved but are only constructed conceptually 
in the matrix splitting of the conservative elemental interval equations. 
The split pieces of the latter, for which the characteristic equations 
may be regarded as eigenfunctions, are the basis of the nodal point finite 
difference equations that are solved in the conservative variables. Fol- 
lowing the approach of Roe, stable nodal point difference equations are 
constructed by choosing split elemental interval equation pieces that 
provide upwind differencings according to the eigenvalues of the associated 
characteristic equations. 

For application to direct implicit matrix methods it is necessary, 
in the delta form, to at least express the time difference of the spatial 
difference of the fluxes in terms of the spatial difference of the time 
differences of the conservative variables. In flux splitting as the result 
of the homogeneous property the flux can be expressed as the product of 
the Jacobian matrix and the conservative variable vector. The delta form 
follows. In the present work a further matrix transform has been found 
that exactly expresses the elemental interval primitive difference equa- 
tions in terms of differences of the conservative variables. Thus with 
the conservative transform of the primitive equations, it has been pos- 
sible in the present method to express the flux difference as a matrix 
transform of the primitive variable differences, hence the nomer flux 
difference splitting and the delta form also follows. 

Finally the split elemental interval equation pieces corresponding 
to stable upwind characteristic equations to the boundaries are used 
in conjunction with auxiliary boundary conditions according to the 
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approach of Kentzer to form stable boundary point difference relations 
that are wholly consistent with the interior point conservative formulation. 

Results of two model problems for which computations have previously 
been reported in the literature have been presented. The solution for 
the first problem, the Burgers equation ramp expansion, found effectively 
the exact solution with no evidence of oscillation or disconnectedness 
(expansion shocks) developing. The solution for the second model problem, 
Shubin's nozzle with supersonic inflow and subsonic outflow has also been 
found to be effectively exact including the shock location. 

VI. Conclusions 

A stable conservative finite difference method has been constructed 
that supports approximately the correct wave propagation and domain of 
dependence of characteristic methods. The splitting is sufficiently 
general to lend itself to a variety of numerical schemes, explicit or 
implicit, iterative or direct, for marching in time or space for the 
compressible Euler and Navier-Stokes equations. Based on preliminary 
results of model problems the upwind split method has the potential for 
both improved robustness and accuracy relative to the present generation 
of unsplit methods. 
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Appendix A 


Supra-Characteristic Matrices 

We couch the discussion in terms of the 2-D axisymmetric 3-D flow 
equations which embody as elements the essential mathematical richness 
of all other cases of interest, i.e. 1-D, quasi 1-D, 2-D and 3-D. For 
a fixed computational mesh based on general curvilinear coordinates, 
the discrete finite volume difference analog of the maximally, conserv- 
ative partial differential equation can be written 

J 6^q + 6^(pW^) + 5^(pW’^) = 0 


J 6^pu + 6^(p u + i^p) + u W'^ + n^P) = 0 

j 6.pV + 6 (p V W^) + i 5 p + 6 (p V W’^) + n 6 p = 0 

t s y s y n 

J 6^E + 6^((E + p)W? + 6^((E + p)W’^) = 0 

where E = + lp(u^ + v^) . 


(Al) 


In the equations = i u + i V , and W’^ = n u + n v , are the cell 

X y X y 

face area weighted, contravariant velocity components, and the metric 
quantities ~ ^ ^ \ ^ ~ \ ^ ^ 

are the x and y projections of the mesh cell face areas, and 
J = 6^^ 'n'Z \ volume of the associated computational 

cell. We note here the v momentum equation does not admit a strong 


Al 



conservation law for the pressure; and we see no reason to introduce 
a weak one with the associated source term that may compromise numerical 
stability. 

It is interesting that the quasi 1-D stream tube equations also 
have the same property of failing to admit a strong conservation law 
for the pressure term of the momentum equation. Indeed a valid repre- 
sentation of the latter set of equations can be obtained from the set 
A1 by dropping the n coordinate differenced terms and the u terms 
including the entire u momentum equation, then replacing v by u 
and by in the remainder, of course the pure 1-D, 2-D, and 
3-D flow formulations admit strong conservation laws like the u momen- 
tum equation of Al, for all the momentum equations of their respective 
sets. 

The gasdynamic equations Al are valid for equilibrium and nonequi- 
librium real gases as well as perfect gases so long as the times for 
relaxing internal atomic and molecular degrees of freedom are short com- 
pared with local residence and reaction times and a local thermodynamic 
equilibrium may be considered to exist. Such real gases (Reference 28) 
are subject to a gas law 


P = p R T (A2) 

derived for a mixture of ideal gases of the real gas composition , 
and the appropriate y for such a system is the specific enthalpy to 
internal energy ratio — . In order to preserve the real gas capability 
of the conservative gasdynamic equations Al in the supra-characteristic 
splittings of derived primitive equations, it is most natural to express 


A2 



; A1 and the derived equations in the primitive variable P = rather 

r 

than the pressure p . 

> Now with the realization fostered in the body of the report that 

we can maintain global conservation while treating each curvilinear 
coordinate direction autonomously, we will restrict the remainder of 
the discussion to splitting the flux difference terms for the- n coordi- 
nate direction. Then the corresponding relations for the c coordinate 
direction can be obtained trivially by replacing the n super or subscripts 
by the appropriate c assuming the associated metric relations. 

With the introduction of the primitive variable P defined above, 
it can be verified by straightforward algebraic manipulation that under 
the matrix multiplication M a vector of primitive equation elemental 
interval differences A^f can be found to identically satisfy the elemental 

interval flux differences of the conservative equation a F . The trans- 

n 

forming matrix M in question is 


1 0 0 0 
u 10 0 

(A3) 

V 0 10 



The right multiplying vector of primitive equation elemental interval 

differences A f is 
n 


•TV 


A3 



7 A A p 

ri n 


pW\u + A (SJy-I)P) 


pw^ V 
n 


n»v 

+ ^ A ((y-1)P) 

V ^ 


(A4) 


Y P n^A^u + Y P a^y + u A^(n/) + vy^P 

and the corresponding vector of elemental interval maximally conservative 

flux differences a F is 

n 

A p W" 
n 

A p u + A (ti„(y"1)P) 

M M ^ 

^ (A5) 

A p V + -^r- A ((y“ 1)P) 
n V n 

A (y P + Hu^ + v2)p 
n 

In the relation and as described in the body of the report, the bar over 
quantities represents the simple average of the values of the quantities 
at either end of the elemental difference interval, (K , K + 1) for n . 
The primitive variable time differences associated with A4 are given by 

•k 

the (column) vector J (6^p , p 6^u , p 6^v , 6^P ) (A6) . These multiplied 
by A3 can be seen to provide measures of the conservative variable time 
differences. 


A4 



Elemental difference supra-characteristic equations can be formed 
of A4 and A6 by multiplying by the matrix T"^ which is 


1 

■p 


yp 


— «L 

pc 

pc 

0 


Hi. 

JL 

pc 

PC 

YP 

1 

|x V 


X 

pc 

PC 

YP 


(A7) 


Here the primed averaged metric quantities are the averaged quantities 
divided by the square root of the sum of their squares or the averaged 
cell face area. In fact, see section III, the primed quantities as 
used in columns 2 and 3 above provide respectively the x and y direc- 
tion cosines of the tangent (row 2) and normal (rows 3 and 4) to the local 
mean n constant coordinate surface. Based on the natural normaliza- 
tions of the primitive equations, the transformation A7 linearizes the 
quasi 1-D elemental difference equations into near orthogonal combinations 
of variables that approximate with first order truncation error scalar 
difference equations. It is well to emphasize here that the intent is 
not to solve such equations in characteristic variables but only to use 
the resulting near orthogonality as an intermediate step to obtain a well 
conditioned matrix formulation in the conservative variables. By row, 
the intermediate equations resulting from A7 support the quasi one dimen- 
sional propagation of local variations in entropy, tangential velocity, 
and the and P~ pressure-normal velocity compatibility relations. 


A5 



r 


The eigenvalues associated with the approximately scalar transformed 

A ^ A A A 

equations are by row + nc , - nc . In the method 

(section II) these eigenvalues remain embedded in the spatial difference 
terms of the primitive equations which are treated exactly. Only the 
sign of the eigenvalues is tested to determine where (forward or backward) 
split pieces of such terms should be "sent", in the sense of Roe, to 
form stable nodal difference equations. 

The inverse T of the matrix T’^ is 


- P 


0 


£ 

2 


£ 

2 


0 


0 



££ ;■ 


- ££' 





(A8) 


0 


0 


YP 

2 


YP 

2 




the primitive variable spatial difference terms be expressed in terms 
of spatial differences of the conservative variables. This end is accom- 
plished exactly by the matrix transform M(Aj^)q = A^f with M(a^) 
given in Table Al. Then with the identity A^F = MA^f , the splittable 
linearizing formulation sketched in relations (1) and (8) of the body 
of the report follows by 


A F = M M“^(a )q = M T D T“^M’^(a )q 
n n n 

where recall D = + D" = I , the identity matrix. 


(A9) 


A6 



Table Al. Conservative variable transform for elemental primitive equation spatial difference terms 


\ - Fvx(-i> 


(Y-l)pn^ _2 

A 

p n 




I 

I 






pu" 


-(Y-l)pn^u 


I -(Y-l)pn 


V Lr 


I 

_L 


I Vx(Y-l) 

I 

_L 


-pW’^ S^^'yV . 

4=r- A - •« A (y-1) 
P n 2 y n 




P V 
-yPh 


— I 

I 

p V n I 


T 


PW»1 1 


7 *1 1 

^ w 

— 1 
-(vrl)p 1 

p V n 1 

A 

V 

lE . - 

p n y 1 

1 


1 

x=3> 

< 

<1 

> 


1 

\V A, 

-H/P 1 

E U V A 1 

P 1 

J 


X - yP 


-=— UA - n. V A 

p n p y n 


-2 72 


V ■ ? \ 


5 ■ “ I W 


YP\ 

A 

p n 


"Y ^ 


-TlvP -2 
U 


- A 



Finally the specification of A , defined in the body of the report 
for the unsteady terms (see equations (10) and (11)) by 

A = M T D r , (AlO) 

is completed here with the matrix inverse of M . The inverse is 

1 0 0 0 

-u 1 0 0 

-V 0 

_2 _2 _ 

U + V -u 

- i(u^ + v^) 


1 0 

-V 1 
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Velocity 

3 0.4 


CSCM First Order Upwind 

00 
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X 

Figure la. Burgers equation solution shown in boxes for the ramp 
expansion by present CSCM method. Exact solution in 
triangles connected by solid line. 



Figure lb. Burgers equation solution for the ramp expansion, reference 27. 
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Figure 2a. Shubins nozzle solution shown in boxes, supersonic inflow, 
subsonic outflow by present CSCM method. Exact solution 
points connected by solid line. 
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Figure 2b. Shubins nozzle, supersonic inflow, subsonic outflow, 
solution as shown in reference 27* 
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Figure 2c. Shubins nozzle, supersonic inflow, subsonic outflow, 
solution as shown in reference 25 . 
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